In MethylKit, genome is fractioned into bins and the methylation levels within each bin are compared using logistic regression (or Fisher’s exact test, if there is only one replicate in each group).
Read in input files. ## Negative control
files.nc <- list.files("./input", pattern = "NC", full.names = T)
files.list.nc <- as.list(files.nc)
anno.nc <- read.table("./input/anno_neg_control.txt.gz", header = T, sep = "\t", stringsAsFactors = F)
myobj.nc <- methRead(
location = files.list.nc, sample.id = as.list(anno.nc$names),
assembly = "hg19", pipeline = "bismarkCoverage", header = F, skip = 0,
dbtype = "tabix", treatment = c(0,1,0,1,0,1), dbdir = "./output/methylDB"
)## creating directory ./output/methylDB ...## Taking input= as a system command ('gunzip -c ./input/NC_rep1.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/NC_rep2.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/NC_rep3.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/NC_rep4.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/NC_rep5.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/NC_rep6.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## compressing the file with bgzip...
## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylBase_4da2c3adde8.txt.bgzfiles.sim <- list.files("./input", pattern = "sim_r", full.names = T)
files.list.sim <- as.list(files.sim)
anno.sim <- read.table("./input/anno_sim_data.txt.gz", header = T, sep = "\t", stringsAsFactors = F)
myobj.sim <- methRead(
location = files.list.sim, sample.id = as.list(anno.sim$names),
assembly = "hg19", pipeline = "bismarkCoverage", header = F, skip = 0,
dbtype = "tabix", treatment = c(0,1,0,1,0,1), dbdir = "./output/methylDB"
)## Taking input= as a system command ('gunzip -c ./input/sim_rep1.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/sim_rep2.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/sim_rep3.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/sim_rep4.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/sim_rep5.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## Taking input= as a system command ('gunzip -c ./input/sim_rep6.cov.gz') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.## compressing the file with bgzip...## making tabix index...## compressing the file with bgzip...
## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylBase_4da22fc53887.txt.bgz## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep1_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep2_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep3_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep4_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep5_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/NC_rep6_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## compressing the file with bgzip...
## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylBase_4da22d8aa5cf.txt.bgz## two groups detected:
## will calculate methylation difference as the difference of
## treatment (group: 1) - control (group: 0)## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylDiff_4da22d8aa5cf.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylDiff_4da22d8aa5cf.txt.bgz.bgzmyDiff.GR.nc <- as(myDiff.nc,"GRanges")
dmr.nc <- data.frame(myDiff.GR.nc)
dmr.nc <- dmr.nc[order(dmr.nc$qvalue),]
colnames(dmr.nc)[7] <- "qval"
# filter regions with mean methylation difference at least 10 % and at least 10 CpGs long
dmr.nc <- dmr.nc[abs(dmr.nc$meth.diff)>=10 && dmr.nc$width >= 10,]
write.table(dmr.nc, file = gzfile("./output/methylKit_dmr_negative_control.txt.gz"),
sep = "\t", row.names = F, quote = F)## checking wether tabix file already exists:## ./output/methylDB/sim_rep1_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## ./output/methylDB/sim_rep2_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## ./output/methylDB/sim_rep3_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## ./output/methylDB/sim_rep4_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## ./output/methylDB/sim_rep5_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## checking wether tabix file already exists:## ./output/methylDB/sim_rep6_tiled.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## compressing the file with bgzip...
## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylBase_4da272d282d3.txt.bgz## two groups detected:
## will calculate methylation difference as the difference of
## treatment (group: 1) - control (group: 0)## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylDiff_4da272d282d3.txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...## flatfile located at: /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylDiff_4da272d282d3.txt.bgz.bgz# filter regions with mean methylation difference at least 10 %
myDiff.sim <- getMethylDiff(myDiff.sim, difference=10)## checking wether tabix file already exists:## /mnt/IM/DKT/courses/sta426-project-dmr-comparison/exp/20181125-differential_analysis_regions-02/20181206-06_methylKit/output/methylDB/methylDiff_4da272d282d3.bgz _type .txt.bgz## tabix file is new.## continuing now ...## compressing the file with bgzip...## making tabix index...myDiff.GR.sim <- as(myDiff.sim,"GRanges")
dmr.sim <- data.frame(myDiff.GR.sim)
dmr.sim <- dmr.sim[order(dmr.sim$qvalue),]
colnames(dmr.sim)[7] <- "qval"
# filter regions with mean methylation difference at least 10 % and at least 10 CpGs long
dmr.sim <- dmr.sim[abs(dmr.sim$meth.diff)>=10 && dmr.sim$width >= 10,]
write.table(dmr.sim, file = gzfile("./output/methylKit_dmr_sim_data.txt.gz"),
sep = "\t", row.names = F, quote = F)Identified DMRs with q-value not higher than 0.05. ## Negative control
c1 <- dmr.nc$qval[dmr.nc$qval <= 0.05]
datatable(
dmr.nc,
rownames = F,
filter = "top", extensions = c("Buttons", "ColReorder"), options = list(
pageLength = 10,
buttons = c("copy", "csv", "excel", "pdf", "print"),
colReorder = list(realtime = FALSE),
dom = "fltBip"
)
) %>%
formatStyle(
"qval",
backgroundColor = styleEqual(c1, rep("#66C2A5", length(c1)))
)## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.htmlc2 <- dmr.sim$qval[dmr.sim$qval <= 0.05]
datatable(
dmr.sim,
rownames = F,
filter = "top", extensions = c("Buttons", "ColReorder"), options = list(
pageLength = 10,
buttons = c("copy", "csv", "excel", "pdf", "print"),
colReorder = list(realtime = FALSE),
dom = "fltBip"
)
) %>%
formatStyle(
"qval",
backgroundColor = styleEqual(c2, rep("#66C2A5", length(c2)))
)## ─ Session info ──────────────────────────────────────────────────────────
## setting value
## version R version 3.5.1 (2018-07-02)
## os Ubuntu 16.04.5 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Zurich
## date 2019-01-10
##
## ─ Packages ──────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.1)
## backports 1.1.3 2018-12-14 [1] CRAN (R 3.5.1)
## bbmle 1.0.20 2017-10-30 [1] CRAN (R 3.5.1)
## bindr 0.1.1 2018-03-13 [1] CRAN (R 3.5.1)
## bindrcpp 0.2.2 2018-03-29 [1] CRAN (R 3.5.1)
## Biobase 2.40.0 2018-10-04 [1] Bioconductor
## BiocGenerics * 0.26.0 2018-10-04 [1] Bioconductor
## BiocParallel 1.14.2 2018-07-08 [1] Bioconductor
## Biostrings 2.48.0 2018-10-29 [1] Bioconductor
## bitops 1.0-6 2013-08-17 [1] CRAN (R 3.5.1)
## callr 3.1.1 2018-12-21 [1] CRAN (R 3.5.1)
## cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.1)
## coda 0.19-2 2018-10-08 [1] CRAN (R 3.5.1)
## colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.1)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
## crosstalk 1.0.0 2016-12-21 [1] CRAN (R 3.5.1)
## data.table * 1.11.8 2018-09-30 [1] CRAN (R 3.5.1)
## DelayedArray 0.6.6 2018-09-11 [1] Bioconductor
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.1)
## devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.1)
## digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.1)
## dplyr 0.7.8 2018-11-10 [1] CRAN (R 3.5.1)
## DT * 0.5 2018-11-05 [1] CRAN (R 3.5.1)
## emdbook 1.3.10 2018-05-19 [1] CRAN (R 3.5.1)
## evaluate 0.12 2018-10-09 [1] CRAN (R 3.5.1)
## fastseg 1.26.0 2018-10-04 [1] Bioconductor
## fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.1)
## GenomeInfoDb * 1.16.0 2018-10-04 [1] Bioconductor
## GenomeInfoDbData 1.1.0 2018-10-04 [1] Bioconductor
## GenomicAlignments 1.16.0 2018-10-04 [1] Bioconductor
## GenomicRanges * 1.32.7 2018-09-20 [1] Bioconductor
## ggplot2 3.1.0 2018-10-25 [1] CRAN (R 3.5.1)
## glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.1)
## gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.1)
## gtools 3.8.1 2018-06-26 [1] CRAN (R 3.5.1)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.1)
## htmlwidgets 1.3 2018-09-30 [1] CRAN (R 3.5.1)
## httpuv 1.4.5.1 2018-12-18 [1] CRAN (R 3.5.1)
## IRanges * 2.14.12 2018-09-20 [1] Bioconductor
## jsonlite 1.6 2018-12-07 [1] CRAN (R 3.5.1)
## knitr 1.21 2018-12-10 [1] CRAN (R 3.5.1)
## later 0.7.5 2018-09-18 [1] CRAN (R 3.5.1)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.5.1)
## lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.1)
## limma 3.36.5 2018-09-20 [1] Bioconductor
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
## MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.1)
## Matrix 1.2-15 2018-11-01 [1] CRAN (R 3.5.1)
## matrixStats 0.54.0 2018-07-23 [1] CRAN (R 3.5.1)
## mclust 5.4.2 2018-11-17 [1] CRAN (R 3.5.1)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.1)
## methylKit * 1.6.3 2018-10-12 [1] Bioconductor
## mime 0.6 2018-10-05 [1] CRAN (R 3.5.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
## numDeriv 2016.8-1 2016-08-27 [1] CRAN (R 3.5.1)
## pillar 1.3.1 2018-12-15 [1] CRAN (R 3.5.1)
## pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.1)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.1)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1)
## processx 3.2.1 2018-12-05 [1] CRAN (R 3.5.1)
## promises 1.0.1 2018-04-13 [1] CRAN (R 3.5.1)
## ps 1.3.0 2018-12-21 [1] CRAN (R 3.5.1)
## purrr 0.2.5 2018-05-29 [1] CRAN (R 3.5.1)
## qvalue 2.12.0 2018-11-28 [1] Bioconductor
## R.methodsS3 1.7.1 2016-02-16 [1] CRAN (R 3.5.1)
## R.oo 1.22.0 2018-04-22 [1] CRAN (R 3.5.1)
## R.utils 2.7.0 2018-08-27 [1] CRAN (R 3.5.1)
## R6 2.3.0 2018-10-04 [1] CRAN (R 3.5.1)
## Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.1)
## RCurl 1.95-4.11 2018-07-15 [1] CRAN (R 3.5.1)
## remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.1)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.1)
## rlang 0.3.0.1 2018-10-25 [1] CRAN (R 3.5.1)
## rmarkdown 1.11 2018-12-08 [1] CRAN (R 3.5.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.1)
## Rsamtools 1.32.3 2018-08-22 [1] Bioconductor
## rtracklayer 1.40.6 2018-08-31 [1] Bioconductor
## S4Vectors * 0.18.3 2018-10-04 [1] Bioconductor
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.1)
## shiny 1.2.0 2018-11-02 [1] CRAN (R 3.5.1)
## stringi 1.2.4 2018-07-20 [1] CRAN (R 3.5.1)
## stringr 1.3.1 2018-05-10 [1] CRAN (R 3.5.1)
## SummarizedExperiment 1.10.1 2018-10-04 [1] Bioconductor
## testthat 2.0.1 2018-10-13 [1] CRAN (R 3.5.1)
## tibble 1.4.2 2018-01-22 [1] CRAN (R 3.5.1)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
## usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.1)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
## xfun 0.4 2018-10-23 [1] CRAN (R 3.5.1)
## XML 3.98-1.16 2018-08-19 [1] CRAN (R 3.5.1)
## xtable 1.8-3 2018-08-29 [1] CRAN (R 3.5.1)
## XVector 0.20.0 2018-10-04 [1] Bioconductor
## yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1)
## zlibbioc 1.26.0 2018-10-04 [1] Bioconductor
##
## [1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.5
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library